/* 
 * File:   cpuFFT.cpp
 * Author: cbender
 * 
 * Created on 5. September 2014, 16:21
 */

#include "cpuFFT.h"


cpuFFT::cpuFFT() {
    gpuInput = (gpu_input_data*) malloc(sizeof(gpu_input_data));
}

cpuFFT::~cpuFFT() {
    free(gpuInput);
}

void cpuFFT::setInput(input_data* input) {
    this->input = input;
}

double cpuFFT::getGreatestAbsNegativeValue() {
    double result = 0;
    for (long i=0; i<this->gpuInput->numPoints; i++) {
        if(frequencyData[i].real<result) {
            result=frequencyData[i].real;
        }
    }
    return fabs(result);
}

void cpuFFT::calculations() {
    __float128 carbonWeight = 0;
    __float128 integerCarbonWeight = 0;
    __float128 carbonVariance = 0;
    __float128 integerCarbonVariance = 0;
    __float128 heliumWeight = 0;
    __float128 integerHeliumWeight = 0;
    __float128 heliumVariance = 0;
    __float128 integerHeliumVariance = 0;
    __float128 nitrogenWeight = 0;
    __float128 integerNitrogenWeight = 0;
    __float128 nitrogenVariance = 0;
    __float128 integerNitrogenVariance = 0;
    __float128 oxygenWeight = 0;
    __float128 integerOxygenWeight = 0;
    __float128 oxygenVariance = 0;
    __float128 integerOxygenVariance = 0;
    __float128 phosphorusWeight = 0;
    __float128 integerPhosphorusWeight = 0;
    __float128 phosphorusVariance = 0;
    __float128 integerPhosphorusVariance = 0;
    __float128 sulfurWeight = 0;
    __float128 integerSulfurWeight = 0;
    __float128 sulfurVariance = 0;
    __float128 integerSulfurVariance = 0;
    __float128 tempMolecularVariance = 0;
    __float128 tempIntegerMolecularVariance = 0;

    // Could be done with loops but should be faster hard-coded

    if (input->numC != 0) {
        carbonWeight = (C12_MASS * input->C12 + C13_MASS * input->C13);
        integerCarbonWeight = (C12_INT_MASS * input->C12 + C13_INT_MASS * input->C13);
        gpuInput->molecularWeight += input->numC * carbonWeight;
        gpuInput->integerMolecularWeight += input->numC * integerCarbonWeight;
        gpuInput->minMolecularWeight += input->numC * C12_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numC * C12_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numC * C13_INT_MASS;
        if (gpuInput->maxIntIsoShiftAsInt < C13_INT_MASS - C12_INT_MASS)
            gpuInput->maxIntIsoShiftAsInt = C13_INT_MASS - C12_INT_MASS;
        carbonVariance += (C12_MASS - carbonWeight)*(C12_MASS - carbonWeight) * input->C12;
        carbonVariance += (C13_MASS - carbonWeight)*(C13_MASS - carbonWeight) * input->C13;
        integerCarbonVariance += (C12_INT_MASS - integerCarbonWeight)*(C12_INT_MASS - integerCarbonWeight) * input->C12;
        integerCarbonVariance += (C13_INT_MASS - integerCarbonWeight)*(C13_INT_MASS - integerCarbonWeight) * input->C13;
    }

    if (input->numH != 0) {
        heliumWeight = (H1_MASS * input->H1 + H2_MASS * input->H2);
        integerHeliumWeight = (H1_INT_MASS * input->H1 + H2_INT_MASS * input->H2);
        gpuInput->molecularWeight += input->numH * heliumWeight;
        gpuInput->integerMolecularWeight += input->numH * integerHeliumWeight;
        gpuInput->minMolecularWeight += input->numH * H1_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numH * H1_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numH * H2_INT_MASS;
        if (gpuInput->maxIntIsoShiftAsInt < H2_INT_MASS - H1_INT_MASS)
            gpuInput->maxIntIsoShiftAsInt = H2_INT_MASS - H1_INT_MASS;
        heliumVariance += (H1_MASS - heliumWeight)*(H1_MASS - heliumWeight) * input->H1;
        heliumVariance += (H2_MASS - heliumWeight)*(H2_MASS - heliumWeight) * input->H2;
        integerHeliumVariance += (H1_INT_MASS - integerHeliumWeight)*(H1_INT_MASS - integerHeliumWeight) * input->H1;
        integerHeliumVariance += (H2_INT_MASS - integerHeliumWeight)*(H2_INT_MASS - integerHeliumWeight) * input->H2;
    }

    if (input->numN != 0) {
        nitrogenWeight = (N14_MASS * input->N14 + N15_MASS * input->N15);
        integerNitrogenWeight = (N14_INT_MASS * input->N14 + N15_INT_MASS * input->N15);
        gpuInput->molecularWeight += input->numN * nitrogenWeight;
        gpuInput->integerMolecularWeight += input->numN * integerNitrogenWeight;
        gpuInput->minMolecularWeight += input->numN * N14_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numN * N14_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numN * N15_INT_MASS;
        if (gpuInput->maxIntIsoShiftAsInt < N15_INT_MASS - N14_INT_MASS)
            gpuInput->maxIntIsoShiftAsInt = N15_INT_MASS - N14_INT_MASS;
        nitrogenVariance += (N14_MASS - nitrogenWeight)*(N14_MASS - nitrogenWeight) * input->N14;
        nitrogenVariance += (N15_MASS - nitrogenWeight)*(N15_MASS - nitrogenWeight) * input->N15;
        integerNitrogenVariance += (N14_INT_MASS - integerNitrogenWeight)*(N14_INT_MASS - integerNitrogenWeight) * input->N14;
        integerNitrogenVariance += (N15_INT_MASS - integerNitrogenWeight)*(N15_INT_MASS - integerNitrogenWeight) * input->N15;
    }

    if (input->numO != 0) {
        oxygenWeight = (O16_MASS * input->O16 + O17_MASS * input->O17 + O18_MASS * input->O18);
        integerOxygenWeight = (O16_INT_MASS * input->O16 + O17_INT_MASS * input->O17 + O18_INT_MASS * input->O18);
        gpuInput->molecularWeight += input->numO * oxygenWeight;
        gpuInput->integerMolecularWeight += input->numO * integerOxygenWeight;
        gpuInput->minMolecularWeight += input->numO * O16_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numO * O16_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numO * O18_INT_MASS;
        if (gpuInput->maxIntIsoShiftAsInt < O18_INT_MASS - O16_INT_MASS)
            gpuInput->maxIntIsoShiftAsInt = O18_INT_MASS - O16_INT_MASS;
        oxygenVariance += (O16_MASS - oxygenWeight)*(O16_MASS - oxygenWeight) * input->O16;
        oxygenVariance += (O17_MASS - oxygenWeight)*(O17_MASS - oxygenWeight) * input->O17;
        oxygenVariance += (O18_MASS - oxygenWeight)*(O18_MASS - oxygenWeight) * input->O18;
        integerOxygenVariance += (O16_INT_MASS - integerOxygenWeight)*(O16_INT_MASS - integerOxygenWeight) * input->O16;
        integerOxygenVariance += (O17_INT_MASS - integerOxygenWeight)*(O17_INT_MASS - integerOxygenWeight) * input->O17;
        integerOxygenVariance += (O18_INT_MASS - integerOxygenWeight)*(O18_INT_MASS - integerOxygenWeight) * input->O18;
    }

    if (input->numP != 0) {
        phosphorusWeight = (P31_MASS * input->P31);
        integerPhosphorusWeight = (P31_INT_MASS * input->P31);
        gpuInput->molecularWeight += input->numP * phosphorusWeight;
        gpuInput->integerMolecularWeight += input->numP * integerPhosphorusWeight;
        gpuInput->minMolecularWeight += input->numP * P31_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numP * P31_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numP * P31_INT_MASS;
    }

    if (input->numS != 0) {
        sulfurWeight = (S32_MASS * input->S32 + S33_MASS * input->S33 +
                S34_MASS * input->S34 + S36_MASS * input->S36);
        integerSulfurWeight = (S32_INT_MASS * input->S32 + S33_INT_MASS * input->S33 +
                S34_INT_MASS * input->S34 + S36_INT_MASS * input->S36);
        gpuInput->molecularWeight += input->numS * sulfurWeight;
        gpuInput->integerMolecularWeight += input->numS * integerSulfurWeight;
        gpuInput->minMolecularWeight += input->numS * S32_MASS;
        gpuInput->minIntegerMolecularWeightAsInt += input->numS * S32_INT_MASS;
        gpuInput->maxIntegerMolecularWeightAsInt += input->numS * S36_INT_MASS;
        if (gpuInput->maxIntIsoShiftAsInt < S36_INT_MASS - S32_INT_MASS)
            gpuInput->maxIntIsoShiftAsInt = S36_INT_MASS - S32_INT_MASS;
        sulfurVariance += (S32_MASS - sulfurWeight)*(S32_MASS - sulfurWeight) * input->S32;
        sulfurVariance += (S33_MASS - sulfurWeight)*(S33_MASS - sulfurWeight) * input->S33;
        sulfurVariance += (S34_MASS - sulfurWeight)*(S34_MASS - sulfurWeight) * input->S34;
        sulfurVariance += (S36_MASS - sulfurWeight)*(S36_MASS - sulfurWeight) * input->S36;
        integerSulfurVariance += (S32_INT_MASS - integerSulfurWeight)*(S32_INT_MASS - integerSulfurWeight) * input->S32;
        integerSulfurVariance += (S33_INT_MASS - integerSulfurWeight)*(S33_INT_MASS - integerSulfurWeight) * input->S33;
        integerSulfurVariance += (S34_INT_MASS - integerSulfurWeight)*(S34_INT_MASS - integerSulfurWeight) * input->S34;
        integerSulfurVariance += (S36_INT_MASS - integerSulfurWeight)*(S36_INT_MASS - integerSulfurWeight) * input->S36;
    }

    // To compensate for the Charge, substract the weight of one electron
    gpuInput->molecularWeight -= ELECTRON_MASS;
    gpuInput->integerMolecularWeight -= ELECTRON_MASS;
    gpuInput->minMolecularWeight -= ELECTRON_MASS;
    gpuInput->integerMolecularWeightAsInt = (int) (gpuInput->integerMolecularWeight + 0.5);

    // Calculate the Molecular Variance
    // the square of the distance of the average mass of an element to each isotope,
    // summed up

    tempMolecularVariance += carbonVariance * input->numC;
    tempMolecularVariance += heliumVariance * input->numH;
    tempMolecularVariance += nitrogenVariance * input->numN;
    tempMolecularVariance += oxygenVariance * input->numO;
    tempMolecularVariance += phosphorusVariance * input->numP;
    tempMolecularVariance += sulfurVariance * input->numS;
    gpuInput->molecularVariance = (cl_real_t) tempMolecularVariance;

    tempIntegerMolecularVariance += integerCarbonVariance * input->numC;
    tempIntegerMolecularVariance += integerHeliumVariance * input->numH;
    tempIntegerMolecularVariance += integerNitrogenVariance * input->numN;
    tempIntegerMolecularVariance += integerOxygenVariance * input->numO;
    tempIntegerMolecularVariance += integerPhosphorusVariance * input->numP;
    tempIntegerMolecularVariance += integerSulfurVariance * input->numS;
    gpuInput->integerMolecularVariance = (cl_real_t) tempIntegerMolecularVariance;
}

void cpuFFT::calculateMassRange() {
    // Calculate the Mass Range
    // This is the smallest Power of two that is still bigger than the
    // molecular Variance wo computed, and it neets to be a power of two for the fft
    cl_ulong dPoints;
    //This is insufficient without adding the one to the end
    dPoints = (cl_ulong) (sqrt(1 + gpuInput->molecularVariance)*10)+0.5;
    gpuInput->bitSize=0;
    while(dPoints){                       // until we found the first bit.
        dPoints >>=1;
        gpuInput->bitSize++;
    }
    gpuInput->massRange=1;
    gpuInput->massRange<<=gpuInput->bitSize;
    // Enhance Resolution. Powers of 2 run well, everything else is not working
    gpuInput->bitSize += gpuInput->resolution;
    gpuInput->numPoints=1;
    gpuInput->numPoints<<=gpuInput->bitSize;
}

void cpuFFT::calculateFrequency() {
    cl_real_t real, imag, freq, theta, r;
    
    // Frequency Domain calculation for every Point in the output array!
    // Normally one would just compute the Mass Range and then do the permutation of
    // data for the fft (reverse the bit order of the input array), but since it will run on the
    // gpu, we will compute the data of the bit reversed array field and put it in the corect one.
    // this will make the overall computation faster since no additional memory access is needed to swap values
    for (cl_ulong i2 = 0; i2 < (cl_ulong) gpuInput->numPoints; i2++) { // has to be unsigned, bitshift will not work as expected otherwise
        // reverse long completely, then shift to the right
        cl_ulong temp = i2;     // temporary variable needed to do the shift
        cl_ulong i =  i2;       // i will be reversed bits of temp
        // reverse order of every bit
        int s = sizeof(i) * CHAR_BIT-1; // extra shift needed at end
        // last bit of i is already correct
        for (temp >>= 1; temp; temp >>= 1)  //as long as temp is not zero
        {   
            i <<= 1;                        // shift bits of i to the left
            i |= temp & 1;                  // if lsb of temp is 1, add it to the right of i
            s--;                            
        }
        i <<= s; // shift when v's highest bits are zero
        i >>= sizeof(i) * CHAR_BIT - gpuInput->bitSize;   // shift it to the right again,
                                                // because the variable can hold more bits than needed
        
        // The upper and lower half of the Points are handled differently
        if (i2 < (cl_ulong) gpuInput->numPoints / 2) {
            freq = (cl_real_t) i2 / gpuInput->massRange;
        } else {
            freq = -(cl_real_t) (gpuInput->numPoints - i2) / gpuInput->massRange;
        }

        //printf("I is %d, I2 is %d\n",i,i2);

        r = 1;
        theta = 0;

        // Do this for every element.

        // Carbon
        if (input->numC) {
            real = input->C12 * cos(2*M_PI * C12_INT_MASS * freq)
                    + input->C13 * cos(2*M_PI * C13_INT_MASS * freq);
            imag = input->C12 * sin(2*M_PI * C12_INT_MASS * freq)
                    + input->C13 * sin(2*M_PI * C13_INT_MASS * freq);
            r *= pow(hypot(real, imag), input->numC);
            theta = theta + input->numC*atan2(imag,real);   //[-pi,pi];	 , taken from midas
        }

        // Helium
        if (input->numH) {
            real = input->H1 * cos(2*M_PI * H1_INT_MASS * freq)
                    + input->H2 * cos(2*M_PI * H2_INT_MASS * freq);
            imag = input->H1 * sin(2*M_PI * H1_INT_MASS * freq)
                    + input->H2 * sin(2*M_PI * H2_INT_MASS * freq);
            r *= pow(hypot(real, imag), input->numH);
              theta = theta + input->numH*atan2(imag,real);   //[-pi,pi];	 , taken from midas
        }

        // Nitrogen
        if (input->numN) {
            real = input->N14 * cos(2*M_PI * N14_INT_MASS * freq)
                    + input->N15 * cos(2*M_PI * N15_INT_MASS * freq);
            imag = input->N14 * sin(2*M_PI * N14_INT_MASS * freq)
                    + input->N15 * sin(2*M_PI * N15_INT_MASS * freq);
            r *= pow(hypot(real, imag), input->numN);
              theta = theta + input->numN*atan2(imag,real);   //[-pi,pi];	 , taken from midas
        }

        // Oxygen
        if (input->numO) {
            real = input->O16 * cos(2*M_PI * O16_INT_MASS * freq)
                    + input->O17 * cos(2*M_PI * O17_INT_MASS * freq)
                    + input->O18 * cos(2*M_PI * O18_INT_MASS * freq);
            imag = input->O16 * sin(2*M_PI * O16_INT_MASS * freq)
                    + input->O17 * sin(2*M_PI * O17_INT_MASS * freq)
                    + input->O18 * sin(2*M_PI * O18_INT_MASS * freq);
            r *= pow(hypot(real, imag), input->numO);
              theta = theta + input->numO*atan2(imag,real);   //[-pi,pi];	 , taken from midas
        }

        // Phosporus
        if (input->numP) {
            real = input->P31 * cos(2*M_PI * P31_INT_MASS * freq);
            imag = input->P31 * sin(2*M_PI * P31_INT_MASS * freq);
            r *= pow(hypot(real, imag), input->numP);
            theta = theta + input->numP*atan2(imag,real);   //[-pi,pi];	 , taken from midas
        }

        // Sulfur
        if (input->numS) {
            real = input->S32 * cos(2*M_PI * S32_INT_MASS * freq)
                    + input->S33 * cos(2*M_PI * S33_INT_MASS * freq)
                    + input->S34 * cos(2*M_PI * S34_INT_MASS * freq)
                    + input->S36 * cos(2*M_PI * S36_INT_MASS * freq);
            imag = input->S32 * sin(2*M_PI * S32_INT_MASS * freq)
                    + input->S33 * sin(2*M_PI * S33_INT_MASS * freq)
                    + input->S34 * sin(2*M_PI * S34_INT_MASS * freq)
                    + input->S36 * sin(2*M_PI * S36_INT_MASS * freq);
            r *= pow(hypot(real, imag), input->numS);
            theta = theta + input->numS*atan2(imag,real);   //[-pi,pi];	 , taken from midas
        }

        /* Convert back to real:imag coordinates and store */
        frequencyData[i].real = r*cos(theta - 2*M_PI*gpuInput->integerMolecularWeightAsInt*freq); //taken from midas
        frequencyData[i].imaginary = r*sin(theta - 2*M_PI*gpuInput->integerMolecularWeightAsInt*freq); //taken from midas
    }

}


void cpuFFT::reset() {
    gpuInput->molecularWeight = 0; // the average molecular weight of the molecule
    gpuInput->minMolecularWeight = 0; // the minimum molecular weight of the molecule, monoMass
    gpuInput->maxIntIsoShiftAsInt = 0; // the maximum difference in Isotope weight as Integer (can be 4,2,1 or 0)
    gpuInput->integerMolecularWeight = 0; // the average molecular weight of the molecule based on the integer values
    gpuInput->integerMolecularWeightAsInt = 0; // the average molecular weight of the molecule based on the integer values
    gpuInput->minIntegerMolecularWeightAsInt = 0; // the minimum molecular weight of the molecule based on the integer values
    gpuInput->maxIntegerMolecularWeightAsInt = 0; // the maximum molecular weight of the molecule based on the integer values
    gpuInput->molecularVariance = 0; // the molecular Variance
    gpuInput->integerMolecularVariance = 0; //the molecular Variance calculated from integer Values
    gpuInput->massRange = 0;
    gpuInput->numPoints = 0;
    gpuInput->resolution = 4;           // multiply the number of points calculated by 2^resolution
    gpuInput->bitSize = 1;
}

void cpuFFT::convertMass() {
    cl_real_t mass, ratio, CorrIntMW;

    if (gpuInput->integerMolecularVariance == 0) ratio = 1;
    else ratio = sqrt(gpuInput->molecularVariance) / sqrt(gpuInput->integerMolecularVariance);
    CorrIntMW = gpuInput->integerMolecularWeight * ratio;

    for (cl_long i = 0; i < gpuInput->numPoints; i++) {
        if (i < gpuInput->numPoints / 2) {
            mass = (double) i / pow(2,gpuInput->resolution) + gpuInput->integerMolecularWeightAsInt;
        } else {
            mass = (double) (i - gpuInput->numPoints) / pow(2,gpuInput->resolution) + gpuInput->integerMolecularWeightAsInt;
        }

        mass *= ratio;
        mass += gpuInput->molecularWeight - CorrIntMW;
        frequencyData[i].imaginary = mass;
    }

}


// only copy the relevant data from the output array back into the list.
// if the weight is too big or too small for the molecule, the result is rejected
// there are also some incorrect peaks intensities from the calculation that should be zero
// and need to be rejected. if a peak is negative, it is obvious that the resulting value is nonsense.
// In order to reject the positive Values that are incorrect, the fact that there are
// as many deviations negative-wise as there are deviations positive-wise is used, and so
// the biggest negative deviation is used as a cut-off value on the positive peak values.
isotope* cpuFFT::copyResult() {
    double cutoffValue = 1.1*getGreatestAbsNegativeValue(); // multiply with 1.1 to get all of them!
    isotope* result = NULL;
    isotope* temp;
    for (cl_long j = 0; j < gpuInput->numPoints; j++) {
        if ((((int) (frequencyData[j].imaginary + 0.5)) >= gpuInput->minIntegerMolecularWeightAsInt) &&
                (((int) (frequencyData[j].imaginary + 0.5)) <= gpuInput->maxIntegerMolecularWeightAsInt) &&
                (frequencyData[j].real > cutoffValue)) {
            temp = (isotope*) malloc(sizeof (isotope));
            temp->mass = frequencyData[j].imaginary;
            temp->p = frequencyData[j].real;
            temp->next = NULL;
            temp->previous = NULL;
            result = add_peak_sorted(result, temp);
            temp = NULL;
        }
    }
    return result;
}

void cpuFFT::fft() {
    cl_int jump;
    cl_int level = 1;
    cl_int i, j, k;
    cl_real_t a, b, c, d, e, f;
    clcomplex swap;

    while ((cl_int) gpuInput->numPoints > level) {
        jump = level << 1;
        a = -1 * (M_PI / level); //inverse FFT
        b = sin(a);
        c = sin(0.5 * a);
        d = -2.0 * c*c;
        e = 1.0;
        f = 0.0;
        for (i = 0; i < level; i++) {
            for (j = i; j < gpuInput->numPoints; j += jump) {
                k = j + level;
                swap.real =         (e * frequencyData[k].real) - (f * frequencyData[k].imaginary);
                swap.imaginary =    (f * frequencyData[k].real) + (e * frequencyData[k].imaginary);
                frequencyData[k].real = frequencyData[j].real - swap.real;
                frequencyData[k].imaginary = frequencyData[j].imaginary - swap.imaginary;
                frequencyData[j].real += swap.real;
                frequencyData[j].imaginary += swap.imaginary;
            }
            c = e;
            e += (c * d) - (f * b);
            f += (f * d) + (c * b);
        }
        level = jump;

    }
}

isotope* cpuFFT::calculate() {
    reset();
    calculations(); // calculate weights and variances
    calculateMassRange(); // calculate the mass Range
    
    frequencyData = (clcomplex*) malloc(gpuInput->numPoints * sizeof(clcomplex));
    
    // This will be done on the CPU
    
    // Take first time measurement
    auto start = std::chrono::high_resolution_clock::now();

    // Calculate the Frequency Data
    calculateFrequency();
    
    // Do the FFT
    fft();
    
    // second time measurement and output of data
    auto finish = std::chrono::high_resolution_clock::now();
    std::cout << std::chrono::duration_cast<std::chrono::nanoseconds>(finish-start).count() << "ns\n";

    // Translate the Mass back to Real Data
    convertMass();
    
    // put the Data into the result
    result = copyResult();

    //Clean up the memory
    free(frequencyData);
    reset(); // set everything back
    return result;
}
